home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / interpolation / interp.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  5.8 KB  |  245 lines
  1. /* interpolation/interp.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Author:  G. Jungman
  21.  */
  22. #include <config.h>
  23. #include <gsl/gsl_errno.h>
  24. #include <gsl/gsl_math.h>
  25. #include <gsl/gsl_interp.h>
  26.  
  27. #define DISCARD_STATUS(s) if ((s) != GSL_SUCCESS) { GSL_ERROR_VAL("interpolation error", (s),  GSL_NAN); }
  28.  
  29. gsl_interp *
  30. gsl_interp_alloc (const gsl_interp_type * T, size_t size)
  31. {
  32.   gsl_interp * interp;
  33.  
  34.   if (size < T->min_size)
  35.     {
  36.       GSL_ERROR_NULL ("insufficient number of points for interpolation type",
  37.                       GSL_EINVAL);
  38.     }
  39.  
  40.   interp = (gsl_interp *) malloc (sizeof(gsl_interp));
  41.   
  42.   if (interp == NULL)
  43.     {
  44.       GSL_ERROR_NULL ("failed to allocate space for interp struct", 
  45.                       GSL_ENOMEM);
  46.     }
  47.   
  48.   interp->type = T;
  49.   interp->size = size;
  50.  
  51.   if (interp->type->alloc == NULL)
  52.     {
  53.       interp->state = NULL;
  54.       return interp;
  55.     }
  56.  
  57.   interp->state = interp->type->alloc(size);
  58.   
  59.   if (interp->state == NULL)
  60.     {
  61.       free (interp);          
  62.       GSL_ERROR_NULL ("failed to allocate space for interp state", GSL_ENOMEM);
  63.     };
  64.     
  65.   return interp;
  66. }
  67.  
  68. int
  69. gsl_interp_init (gsl_interp * interp, const double x_array[], const double y_array[], size_t size)
  70. {
  71.   if (size != interp->size)
  72.     {
  73.       GSL_ERROR ("data must match size of interpolation object", GSL_EINVAL);
  74.     }
  75.   
  76.   interp->xmin = x_array[0];
  77.   interp->xmax = x_array[size - 1];
  78.  
  79.   {
  80.     int status = interp->type->init(interp->state, x_array, y_array, size);
  81.     return status;
  82.   }
  83. }
  84.  
  85. const char *
  86. gsl_interp_name(const gsl_interp * interp)
  87. {
  88.   return interp->type->name;
  89. }
  90.  
  91. unsigned int
  92. gsl_interp_min_size(const gsl_interp * interp)
  93. {
  94.   return interp->type->min_size;
  95. }
  96.  
  97. void
  98. gsl_interp_free (gsl_interp * interp)
  99. {
  100.   if (interp->type->free)
  101.     interp->type->free (interp->state);
  102.   free (interp);
  103. }
  104.  
  105.  
  106.  
  107. int
  108. gsl_interp_eval_e (const gsl_interp * interp,
  109.                    const double xa[], const double ya[], double x,
  110.                    gsl_interp_accel * a, double *y)
  111. {
  112.   if (x < interp->xmin)
  113.     {
  114.       *y = ya[0];
  115.       return GSL_EDOM;
  116.     }
  117.   else if (x > interp->xmax)
  118.     {
  119.       *y = ya[interp->size - 1];
  120.       return GSL_EDOM;
  121.     }
  122.  
  123.   return interp->type->eval (interp->state, xa, ya, interp->size, x, a, y);
  124. }
  125.  
  126. double
  127. gsl_interp_eval (const gsl_interp * interp,
  128.          const double xa[], const double ya[], double x,
  129.          gsl_interp_accel * a)
  130. {
  131.   double y;
  132.   int status = interp->type->eval (interp->state, xa, ya, interp->size, x, a, &y);
  133.  
  134.   DISCARD_STATUS(status);
  135.  
  136.   return y;
  137. }
  138.  
  139.  
  140. int
  141. gsl_interp_eval_deriv_e (const gsl_interp * interp,
  142.                          const double xa[], const double ya[], double x,
  143.                          gsl_interp_accel * a,
  144.                          double *dydx)
  145. {
  146.   if (x < interp->xmin)
  147.     {
  148.       *dydx = 0.0;
  149.       return GSL_EDOM;
  150.     }
  151.   else if (x > interp->xmax)
  152.     {
  153.       *dydx = 0.0;
  154.       return GSL_EDOM;
  155.     }
  156.  
  157.   return interp->type->eval_deriv (interp->state, xa, ya, interp->size, x, a, dydx);
  158. }
  159.  
  160. double
  161. gsl_interp_eval_deriv (const gsl_interp * interp,
  162.                const double xa[], const double ya[], double x,
  163.                gsl_interp_accel * a)
  164. {
  165.   double dydx;
  166.   int status = interp->type->eval_deriv (interp->state, xa, ya, interp->size, x, a, &dydx);
  167.  
  168.   DISCARD_STATUS(status);
  169.  
  170.   return dydx;
  171. }
  172.  
  173.  
  174. int
  175. gsl_interp_eval_deriv2_e (const gsl_interp * interp,
  176.                           const double xa[], const double ya[], double x,
  177.                           gsl_interp_accel * a,
  178.                           double * d2)
  179. {
  180.   if (x < interp->xmin)
  181.     {
  182.       *d2 = 0.0;
  183.       return GSL_EDOM;
  184.     }
  185.   else if (x > interp->xmax)
  186.     {
  187.       *d2 = 0.0;
  188.       return GSL_EDOM;
  189.     }
  190.  
  191.   return interp->type->eval_deriv2 (interp->state, xa, ya, interp->size, x, a, d2);
  192. }
  193.  
  194. double
  195. gsl_interp_eval_deriv2 (const gsl_interp * interp,
  196.                 const double xa[], const double ya[], double x,
  197.                 gsl_interp_accel * a)
  198. {
  199.   double d2;
  200.   int status = interp->type->eval_deriv2 (interp->state, xa, ya, interp->size, x, a, &d2);
  201.  
  202.   DISCARD_STATUS(status);
  203.  
  204.   return d2;
  205. }
  206.  
  207.  
  208. int
  209. gsl_interp_eval_integ_e (const gsl_interp * interp,
  210.                          const double xa[], const double ya[],
  211.                          double a, double b,
  212.                          gsl_interp_accel * acc,
  213.                          double * result)
  214. {
  215.   if (a > b || a < interp->xmin || b > interp->xmax)
  216.     {
  217.       *result = 0.0;
  218.       return GSL_EDOM;
  219.     }
  220.   else if(a == b)
  221.     {
  222.       *result = 0.0;
  223.       return GSL_SUCCESS;
  224.     }
  225.  
  226.   return interp->type->eval_integ (interp->state, xa, ya, interp->size, acc, a, b, result);
  227. }
  228.  
  229.  
  230. double
  231. gsl_interp_eval_integ (const gsl_interp * interp,
  232.                const double xa[], const double ya[],
  233.                        double a, double b,
  234.                gsl_interp_accel * acc)
  235. {
  236.   double result;
  237.   int status = interp->type->eval_integ (interp->state, xa, ya, interp->size, acc, a, b, &result);
  238.  
  239.   DISCARD_STATUS(status);
  240.  
  241.   return result;
  242. }
  243.  
  244.  
  245.